Dynamics of Phase Transitions by Hysteresis Methods I 
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In studies of the QCD deconfining phase transition or crossover by means of heavy ion experi- 
ments, one ought to be concerned about non-equilibrium effects due to heating and cooling of the 
system. Motivated by this, we look at hysteresis methods to study the dynamics of phase transi- 
tions. Our systems are temperature driven through the phase transition using updating procedures 
in the Glauber universality class. Hysteresis calculations are presented for a number of observables, 
including the (internal) energy, properties of Fortuin-Kasteleyn clusters and structure functions. 
We test the methods for 2d Potts models, which provide a rich collection of phase transitions with 
a number of rigorously known properties. Comparing with equilibrium configurations we find a 
scenario where the dynamics of the transition leads to a spinodal decomposition which dominates 
the statistical properties of the configurations. One may expect an enhancement of low energy gluon 
production due to spinodal decomposition of the Polyakov loops, if such a scenario is realized by 
nature. 

PACS numbers: PACS: 05.50.+q, 11.15.Ha, 25.75.-q, 25.75.Nq 



I. INTRODUCTION 

Quantum chromodynamics has well established phase 
transitions in certain limiting cases. In the limit of van- 
ishing quark masses it has the chiral phase transition 
from the phase of broken chiral symmetry at low tem- 
peratures to the chiral symmetric phase at high temper- 
atures. In the limit of infinite quark masses one finds the 
deconfinement transition from the Z(3)-symmetric low 
temperature phase with confinement to the i?(3V_broken 
phase at high temperatures, for a review see For 
physical quark masses m u ,md of the order of 10 MeV 
and m s of the order of 150 MeV it is suggested by lat- 
tice simulations and effective models [J that neither a 
chiral nor a deconfinement transition occurs in the sense 
that there are thermodynamic singularities. 

Lattice gauge theory investigations of the finite tem- 
perature phase transitions of QCD have, with some no- 
table exceptions , been limited to studies of their equi- 
librium properties, whereas in nature these transitions 
are governed by a temperature, or otherwise, driven dy- 
namics. Even when a proper phase transition does not 
exist, a question is whether one may expect observ- 
able remnants of the phase conversion because of off- 
equilibrium effects. 

In the early universe the effects of the dynamics are 
most likely negligible, since the cooling process is deter- 
mined by the Hubble expansion of the universe that is 
slow compared to the typical time scales of strong inter- 
actions, which are of the order of 10 -23 sec. In heavy 
ion collisions this is different. A rapid heating (quench) 
of the nuclei at the "little bang" event is followed by a 
slower cooling process. The lifetime of the emerging sys- 
tem appears to be sufficiently long to equilibrate a phe- 



nomenological quark-gluon plasma 0, |(J, although the 
dynamics in the time period of the phase conversion may 
proceed out of equilibrium. Finite size corrections may 
play a role, because the system is not large compared to 
the typical spatial scale of strong interactions, i.e. 1 fm. 
One should also address the question, whether the initial 
quench could lead to domains of distinct average Z 3 3- 
ality, with interfaces between them, which have relatively 
long relaxation times. 

On an effective level (in the framework of the 0(4)- 
model) one has studied dynamical effects on the chiral 
phase transition 0, Although the largest equilib- 
rium correlation length (that of the pion with a mass of 
w 137 MeV) is not large compared to the intrinsic QCD 
scale (e.g., set by Amom), as result of a quenched cooling 
process one may get disoriented chiral condensates via 
spinodal decomposition. We are interested in the analo- 
gous question for the deconfinement transition 9]. One 
could get a disoriented condensate of Polyakov [lfj loops 
and an associated production of low-momentum gluons. 

Polyakov loops behave effectively like 3c? spin vari- 
ables 0, 0, Q and the Potts-model in three di- 
mensions with q = 3 states gives an effective description 
of the deconfinement transition (more sophisticated spin 
models are also considered 13]). By adding an exter- 
nal field 0| , one can represent the effect of finite quark 
masses. Even this simplification is not yet a suitable ba- 
sis for a numerical investigation. To get confidence in our 
computational methods, we simulated g-state Potts mod- 
els in 2c?, for which a number of rigorous results [prl \v\ 
allow for cross checks. We set the external field to zero 
and choose q = 2, 4, 5 and 10, corresponding to a weak 
second order, a strong second order, a weak first order 
and a strong first order phase transition, respectively. 
The difference between weak and strong second order 
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transitions is explained in section II I II For a review of 
Potts models see Wu |Tg| . 

We use hysteresis methods to investigate the phase 
transition in the Glauber |l9j ] dynamics. The universality 
class of Glauber dynamics, model A in the classification 
of Ref . pp| . contains local Monte Carlo (MC) updating 
schemes which imitate the thermal fluctuations of na- 
ture. Studying the computer time evolution of Glauber 
dynamics gives an overview of a scenario which allows for 
a variation of the speed of the phase transition. Notably, 
the notion of the Minkowskian time is lost in the conven- 
tional quantum field-theoretical formulation of an equi- 
librium ensemble |2l| which is used in numerical simula- 
tions. To study the time evolution of this field-theoretic 
ensemble, one has to find a way to re-introduce a proper 
dynamics. The hope is that the thus generated configu- 
rations are typical for the dynamical process. 

Our observables are the internal energy, properties 
of Fortuin-Kasteleyn (FK) clusters [2jJ , and structure 
functions. The results from equilibrium configurations 
are compared with those from configurations that are 
dynamically driven through the hysteresis cycles. In 
all cases we find that the dynamics induces remarkably 
strong signals for a spinodal decomposition. With in- 
creasing q similar signals become very weak for the equi- 
librium phase transition. 

In the next section we discuss in more detail the basic 
concepts used in this paper. Our numerical investigations 
are reported in section IIIII where subsections deal with 
bulk properties, FK clusters and structure functions. A 
brief summary and conclusions are given in the final sec- 
tion H2] Article II of this series will be devoted to a 
study of the 3d 3-state Potts model in an external mag- 
netic field. 



II. PRELIMINARIES 

Our (computer) time-dependent Hamiltonian is 

H(t) = -0(t)E (1) 

where 

E = -2 <5(r(r,t),(T(r',t) ■ (2) 

(f,r') 

Here the sum runs over all nearest neighbor sites r and 
r ' and a takes the values 1, . . . ,q. In this paper we rely 
on symmetric lattices of N = L x L spins. For suit- 
ably chosen values of /3 m - ln and /3 max , we run the system 
at various cooling/heating rates in cycles from /3 m in to 
Anax and back. Hysteresis methods plaved some role in 
the early days of lattice gauge theory |2J| , but have ap- 
parently been abandoned. Possibly, the reason is that 
one does not learn much from a single hysteresis. How- 
ever, averages over large numbers of heating and cooling 
cycles have to our knowledge not been analyzed in the 
literature. By creating a large number of cycles, ensem- 
ble averages of dynamical configurations are obtained at 



selected temperatures T — 1/(3. For each temperature 
away from the endpoint of the cycles two distinct av- 
erages exist, one on the heating and the other on the 
cooling branch of the cycles. 

The spins are updated by an algorithm which is within 
the Glauber class. Examples are single- and multi- 
hit Metropolis, as well as heat-bath updating methods, 
where the lattice sites may be visited randomly or in some 
systematic order. At critical points (i.e. at second order 
transitions) the slowing down of such algorithms is gov- 
erned by universal exponents. A counterexample is the 
Swendsen-Wang 25] algorithm, which updates entire FK 
clusters. Clearly, such an updating does not correspond 
to thermal fluctuations of nature. The purpose of the 
Swendsen-Wang algorithm is to speed up the dynamics of 
second order phase transitions in computer simulations. 

When driving the system through the transition, the 
phase conversion may be dominated by metastable or 
unstable states of matter. If a system is brought into a 
metastable state, it will be unstable against finite, local- 
ized fluctuation. This scenario is called nucleation. It 
may allow the system to reach a metastable equilibrium 
before a large enough fluctuations occurs. If the sys- 
tem is brought into an unstable state, infinitesimal, non- 
localized amplitude fluctuations lead to an immediate on- 
set of the decay of the unstable state. This scenario is 
called spinodal decomposition. It may lead to long-range 
correlations, in a sense similar to those encountered in 
equilibrium close to second order phase transitions. 

The concept of nucleation as well as the spinodal 
were first introduced by Gibbs as early as 1877, where 
the spinodal was defined as a limit for metastability of 
fluid gases. But only in the late 1950s it became ap- 
parent that a phase beyond the spinodal decomposes 
by diffusional clustering mechanism quite different from 
the nucleation and growth mechanism encountered for 
metastable states. In his classic review |26| Cahn includes 
an account of the historical development. The mod- 
ern theory uses effective diffusional differential equations 
(originally an idea of Hillert 1271') to distinguish dynami- 
cal universality classes, see Ref. |20t EH E^| for reviews. A 
sharp distinction between infinitesimal (spinodal) and fi- 
nite (nucleation) fluctuations is, strictly speaking, a mean 
field concept. In real systems, where fluctuations are im- 
portant, the boundary separating nucleation from spin- 
odal decomposition is not perfectly sharp. 

The numerical investigations, we are aware off, inves- 
tigate the spinodal versus the nucleation scenario after a 
quench, which may either lead into the metastable region 
(nucleation) or beyond it (spinodal decomposition). See 
Miller and Ogilvie |j] in the context of lattice gauge the- 
ory. Our hysteresis approach differs in this respect. The 
continued change of the external temperature prevents 
the system from ever reaching equilibrium, but implies on 
the other hand a smoother dynamics, because the tem- 
perature changes only in small steps. Under laboratory 
conditions there is never a perfect quench and in some 
situations our hysteresis approach may allow to model 
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laboratory condition more realistically than a quench. 
We measured many observables in each hysteresis cycle. 
In this paper we report selected results for the energy 
FK clusters and structure functions. In more detail the 
data will be analyzed and presented in Ref. [3(j • 

We measure FK clusters instead of geometrical clus- 
ters, because their statistical definition accounts for the 
fact that neighboring spins may not only be aligned by 
the spontaneous magnetization but also by random fluc- 
tuations. It is only then that the Kertesz j^lj line of per- 
colation coincides with the phase transition, see Ref.|3^| 
for a review of this and related topics. In contrast to the 
stochastic definition, the geometric definition connects 
aligned spins with certainty and leads to an overcounting 
of ordered clusters. While the FK works well for Potts 
models, a generalization to gauge theories is not known. 
This is closely related to the fact that a cluster updating 
algorithm is not known for gauge theories. 

We are interested in the effects of dynamic heating and 
cooling on the cluster structure, in particular in the ques- 
tion, whether one may still find observable signals, even 
when their is no longer a transition in the strict ther- 
modynamic sense. There are similarities and differences 
to the program of Satz [3j|. Satz focuses on geometric 
properties of FK clusters and would like to extract from 
their equilibrium distribution signals for the phase con- 
version when there is no proper phase transition. Wc 
are trying to find signals for the phase conversion due 
to the deviations from equilibrium. For nucleation one 
expects compact clusters, due to the non-zero interfacial 
tension between the ordered and the disordered phase. 
For spinodal decomposition clusters of each of the or- 
dered states will grow unrestricted by such an interfacial 
tension, building domain walls between the distinct or- 
dered states. For nucleation we expect the maximum 
cluster surface to grow to a size cL d ~ x with c ss 2 for 
strong first oder transitions (c = 2 for the smallest surface 
of a cluster which percolates). For spinodal decomposi- 
tion we expect considerably larger values, comparable to 
the largest values one finds on equilibrium configurations 
in the neighborhood of a second order phase transition. 

In our simulations we record the following cluster ob- 
servables: their number, the mean volume, the maximum 
volume, the mean surface area, the maximum surface 
area, the gyration radius and the percolation probability. 
The volume of a cluster is simply the number of spins it 
contains. The cluster surface is defined on the links of 
the dual lattice, which corresponds to the d — 1 dimen- 
sional hyperspace of the original lattice. The percolation 
probability p is the probability to find at least one clus- 
ter that percolates. For our periodic lattices this means 
that the cluster connects to itself through the boundary 
conditions, in any one of the two directions. 

We analyze the structure function in momentum space 
for signals of spinodal decomposition. Let m qo — 
(fia(f,t),q ) denote the magnetization in direction go £ 
{1, ...,q}. By introducing a Potts spin S qo (r,t) — 



b<r(r,t).q we can write the correlation function 



g(r,?,t) = (S a ^ t ).a(r',t)} 
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in the familiar form 



9-1 9-1 
g(r,7,t) = {S qo (r,t)S qo (7,t)) £ {S qo f 

9o=0 90=0 



(3) 



(4) 



The structure factor (function) is the Fourier transform 
of the correlation function 



S(M) = ^X>(|£|,*)exp[tfci2| 



(5) 
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where R = r — f . Some straightforward algebra trans- 
forms this into 



8(k,t) 




^ <W,*),go exp[ifcr| 
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This is simply the time-dependent version of the equi- 
librium structure factor. In condensed matter experi- 
ments the magnitude of the structure function is directly 
observable in X-ray neutron and light scattering exper- 
iments, compare, e.g., Ref. |34j|. Unfortunately, it ap- 
pears to us that direct measurements in high energy ex- 
periments are unrealistic. In our simulations we expect 
pronounced peaks (similar as for equilibrium configura- 
tion near second order phase transitions) for S(k, t) in 
the case of a phase conversion by spinodal decomposi- 
tion and no such signals in the case of a conversion by 
nucleation and growth. 



III. NUMERICAL RESULTS 

The data presented in this paper rely on systematic 
updating for which the Potts spins are updated in se- 
quential order, each spin once during one sweep. We 
did a number of cross-checks using random updating for 
which the spins are updated in random order, in the av- 
erage each spin once during one sweep. Besides a slowing 
down of the dynamics by a factor of about 0.6 for ran- 
dom updating, we observed no noticeable changes of the 
results checked. 

The temperature (3 — 1/T is changed by ±A/3 af- 
ter every sweep (we experimented also with temperature 
changes after each spin update and found no differences 
within our statistical errors). Our stepsize A/3 is propor- 
tional to the inverse volume of the system 



A/3 = 



2 (/3max ftmin ) 

np L 2 



(7) 
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where /3 m i n and /3 max define the terminal temperatures 
and the integer np = 1,2,... is varied. Equilibrium con- 
figurations are recovered in the limit nn — > oo (A/3 = 0). 
In nature the fluctuations per spin per time unit (here 
the unit of one MC sweep) set the scale for the dynam- 
ics. Our choice of A/3 is motivated by our interest in the 
question whether a dynamics, which slows down with 
volume size may still dominate the nature of the tran- 
sition. Relying on the heat-bath method, each of our 
systems is driven through at least 640 cycles, each start- 
ing from an equilibrated, disordered configuration. Er- 
ror bars are calculated with respect to 32 jackknife bins. 
In an exploratory simulation 35] of 2c? Potts models the 
Metropolis algorithm was employed, but it turns out that 
the heat-bath method saves CPU time. 

In this article we simulate the 2d g-state Potts model 
for q = 2, 4, 5 and 10. This allows us to compare the in- 
fluence of the Glauber dynamics for a weak second order, 
a strong second order, a weak first order and a strong first 
order phase transition. Our terminology "strong second 
order phase transition" may need some explanation. For 
a finite system of volume L d the partition functions is a 
polynomial in u — exp(— /3) that takes positive values on 
the real axis. For first and second order phase transitions 
the imaginary part of the partition function's zero closest 
to the real axis scales like u® ~ L^ 1 / 1 ' , where v = 1/d for 
a first order transition and 1/d < v < 2/d for a second 
order transition. The fluctuations of the energy are gov- 
erned by the exponent a of the specific heat for which we 
assume the hyperscaling relation 36] a = 2 — dv. There- 
fore, a = 1 for first order transitions and < a < 1 for 
second order transitions. To determine the implications 
for the finite size scaling of the energy fluctuations, we 
use the link expectation value of the energy 
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et((3) = (5 a( r, t )Mr',t)) = -{E)/(2dL d ) 



(8) 



where r and r are nearest neighbor sites. The values of 
ei are conveniently located in the range < e; < 1 with 
e;(0) = 1/q and e;(oo) = 1. To leading order in L, finite 
size scaling theory predicts the fluctuation of e; to scale 
like 



({ ei f) - ( ei ) 2 ~ L«/"" 



(9) 



for /3 at the transition point (3 C = 1/T C . For first or- 
der phase transitions a = 1 holds and the left-hand-side 
of equation 10 approaches a finite value, proportional 
to the square of the latent heat Ae;. For second order 
phase transitions the left-hand-side scales to zero. In this 
sense a weak second order transition is one with a close 
to zero or a = and a cusp or logarithmic singularity, 
while a strong second order transition has a close to one. 
First order transitions are weak when Ae; <C 1 holds 
and strong when Ae; becomes of order one, say from 
Aei > 0.1 on. For our choices of q the analytical val- 
ues of /3 C , a and Ae; are compiled in tabled Our 
values of /3 m ; n and /3 max for equation Q and a numerical 
result, Ae;, as explained in the following subsection, are 
also given in this table. 



TABLE I: The (infinite volume) phase transition tempera- 
tures /3 C = 1/Tc, the specific heat exponent a and the latent 
heats of selected g-state Potts models in two dimensions. For 
the latent heats the negative energy per link Ae; is given and 
Ae; is an estimate from hysteresis cycles. 
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FIG. 1: Energy iJHJ hysteresis curves for np 
right: q = 2, 4, 5, and 10. 



1. From left to 



In steps of 20 our lattice sizes range from L — 20 to 
L = 100. For the smaller systems all hysteresis runs are 
done on a single PC, while for the larger lattices up to 
32 PCs are used, dividing our entire run in 32 bins of 
at least 20 hysteresis loops each. In each case a short 
equilibrium run of q L 2 sweeps was initially performed 
at /3 m j n , where the systems equilibrate easily, because 
they are highly disordered. For comparison with equilib- 
rium configurations we performed multicanonical [s?! Issl ] 
(see the next subsection) as well as conventional, canon- 
ical simulations. The reason for the conventional canon- 
ical equilibrium simulations is that one needs to know 
the temperature to generate FK clusters. They were 
performed at many temperatures and in each case 640 
measurements were taken after at least 20g L 2 sweeps for 
reaching equilibrium. 



A. Internal Energy 

For a first order phase transition, the slowing down of 
the canonical equilibrium Markov process is exponential 
in computer time, ~ exp[2 f 8 L d ~ 1 ] 1 where f s is the in- 
terfacial tension (see [l7| for the analytical values). In 
this case we expect an energy hysteresis to survive in the 
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shows the thus obtained estimates Ae;(L) together with 
fits of the form 



20 30 40 50 60 70 80 90 100 
L 

FIG. 2: Latent heat estimates from np = 1 hysteresis curves. 



limit L — > oo and A/3(L; riff) — * for any fixed value of 
rip in equation (JJJ). The shape of the hysteresis can then 
be used to define finite volumes estimators of physical 
variables, such as the transition temperature and the la- 
tent heat. The infinite volume limits of these estimators 
are supposed to be independent of any fixed choice of rip. 

For the second order phase transitions of the q = 2 and 
q = 4 models the analysis is more subtle. The Markov 
process slows only down like L z with z rs 2 There- 
fore, one still expects a hysteresis in the limit L — > oo 
and np fixed, only the opening has no longer the inter- 
pretation of a finite volume estimator of the equilibrium 
latent heat. A finite size scaling analysis of the hystere- 
sis as function of np(L) should allow to identify second 
order transitions. This analysis is not pursued here. 

For rip = 1 and selected lattice sizes we show in figure^ 
our energy ijHJ) hysteresis data. The ordinate is scaled to 



C(q) [e0)-~ 

q 



with C{q) 



(io) 



so that, independently of q, the range [0, 1] gets covered 
when (3 is varied from — ^ oo. If one wants to compare 
the present heat bath with Metropolis results j^H, the 
better efficiency of the heat bath algorithm is such that 



Metropolis 



g^hcatbath Qyg^ fa De use d. From left to 



right in figure ^hysteresis loops for the cases q = 2, 4, 
5 and 10 are visible. For clarity of the figure we have 
omitted error bars and for q = 4 and 5 also the L = 40 
and 60 lattices. Notable is that the hysteresis curves for 
the q — 4 strong second order transition and the q = 5 
weak first order transition are quite similar. From q = 10 
to q = 2 there is a gradual, not an abrupt, deformation 
of the shape of the hysteresis. 

To analyze the physical content of the hysteresis curves 
of figure ^ in more detail, we define the finite volume es- 
timators of the inverse transition temperature (3 C (L) and 
of the latent heat Aei(L) by their values at the maximum 
opening of the corresponding hysteresis curve. Figure |3 



Aez(L) = Aej + j- 



(11) 



where a\ is a constant. For q = 10 the left-hand-side or- 
dinate applies and for the other g-values the right-hand- 
side ordinate. Because of the distinct scales the differ- 
ence between the q = 10 and the q = 5 estimators is 
large, while the general behavior of the fitting curve ap- 
pears to be quite similar for all q-values. The obtained 
infinite volume estimates Aej are given in table HI For 
q = 10 the estimate is in excellent agreement with the 
analytical result, but this is not at all the case for the 
other g-values. Instead, the q = 2, 4 and 5 estimates 
overshoot the equilibrium values considerably. This does 
not come as a surprise, because we already noted that, in 
the infinite volume limit and for fixed rip, a finite open- 
ing of the hysteresis survives even for the second order 
phase transitions. Obviously, the opening has no longer 
the interpretation of an estimator of the equilibrium la- 
tent heat. Instead, the phenomenon illustrates that the 
dynamics tends to wash out differences of the equilibrium 
properties of the transitions. 

Performing a similar analysis for [3 C (L) and compar- 
ing the infinite volume estimates with the analytical re- 
sults, we get accuracies of about ±1% for all q. So, 
we find no problem in locating the equilibrium transi- 
tion temperature from the information of the dynamics. 
The accuracy of these dynamical estimates is not com- 
petitive with the best equilibrium methods. E.g., fitting 
the pseudocritical /3-values of the multicanonical 10-state 
Potts model simulation [37j self-consistently to the form 
P C (L) = P c + c/L 2 gives j3 c = 0.713032 (16) (using our 
energy convention ©). The purpose of our present study 
is not to calculate high precision estimates of equilibrium 
quantities, but to investigate the deviations from equilib- 
rium due to the imposed dynamics. To understand the 
dynamics of our finite volume transitions in more detail, 
we analyze in the next two subsections the behavior of FK 
clusters and structure functions on our configurations. 

In a last remark about the hysteresis curves of the in- 
ternal energy, we like to mention that we have also gener- 
ated equilibrium data for all cases using the multicanon- 
ical method. As expected, the thus obtained e/(/3) func- 
tions fall inside the hysteresis curves of figure ^ Some 
more details are given in Ref . |35| . 



B. Cluster Properties 

We limit our presentation to a few of the cluster ob- 
servables we measure (more details will be given in 30] ). 
The largest cluster surface turns out to be interesting, 
because it exhibits pronounced peaks in the transition 
region. We use the normalization 



S r = 



# surface links of cluster c 
L^ 1 



(12) 
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FIG. 3: The largest cluster surface for the 10-state Potts 
model on various lattice sizes as indicated in the figure (the 
extensions are the value of np and e for equilibrium). 




FIG. 5: The largest cluster surface for the 2-state Potts model 
on various lattice sizes as indicated in the figure. 
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FIG. 4: The probability p of having a percolating cluster for 
the 10-state Potts model on various lattices. 



FIG. 6: The largest cluster surface for equilibrium simula- 
tions. 



for our cluster surfaces. A link which connects a site 
of the cluster with the site of another adjacent cluster 
is defined to be a surface link. The surface links can 
be mapped on the (d— l)-dimensional hypercubes which 
enclose the cluster on the dual lattice. The largest surface 
is simply defined as 

Smax = max-fSc} (13) 

where the maximum is taken over all clusters c of the 
configuration at hand. 

For the 10-state model results for the largest cluster 
surface of the rip — 1 hysteresis cycle are shown in fig- 
ure |3J The arrows indicate the flow of the hysteresis cy- 
cles. During the heating and cooling parts of the cycles, 
the surface areas peak at distinct values, (3 = f3p Cak - This 
is striking evidence that the geometry of the FK clusters 
is distinct during cooling and heating. Due to our use of 
stochastic (in contrast to geometrical) clusters, the equi- 



librium transition temperature value is pinched between 
the temperatures at which the two peaks are located. 

It can be understood that the peaks of S max (f3) are 
related to percolation. For the (3 — > /3 m ax half-cycle the 
picture is that the cluster with the largest surface per- 
colates due to the periodic boundary conditions. Until 
the cluster percolates, its surface area increases, while it 
is decreasing after percolation (as only small islands of 
the false phase remain eventually). Relying on the same 
data as for figure El we show in figure 0] the percolation 
probability p. It is seen that the temperatures of the 
Smax peaks correspond approximately to the steepest in- 
crease/decrease of the percolation probabilities. 

Another observation from figure is that for the half- 
cycle (3 — > /3 m in the peaks of S'max are even more pro- 
nounced than for (3 — > /3 max - This is in accordance with 
a very rapid fall-off of the percolation probability for the 
f3 — » /3 m ; n half-cycle. Our interpretation is that the re- 
sponse to the temperature change is more rapid when the 



7 



system enters the disordered phase than when it enters 
the ordered phase. Such a change in relaxation scales 
may be expected for a strong first order transition (be- 
cause both phases are separated by a gap in the energy 
and not continuously related), while one would expect 
that the response times under heating and cooling are 
similar for a weak second order phase transition. Indeed, 
figure El shows that the two peaks are of almost equal 
height for the 2-state Potts (Ising) model. The S max re- 
sults for the q = 4 and q — 5 models (no figures shown) 
are in-between the two scenarios, but certainly closer to 
q = 10 than to q — 2. The difference between q = 4 and 
q — 5 is minor. 

Also shown in figures El and El are results for S max ([3) 
from equilibrium simulations on L = 80 lattices. They 
are barely visible, because they are to a large extent cov- 
ered by the curves of the j3 — ► /3 m i n half-cycle. Therefore, 
we plot the equilibrium curves for all our g-values sepa- 
rately in figureE! The peaks show a marked increase from 
q = 10 (right) to q — 2 (left). For first order phase tran- 
sitions the interface tension implies that the free energy 
increases with the cluster surfaces. The stronger the first 
order transition is, the more the system tries to minimize 
interfaces. For a second order phase transition there is no 
(disorder-order) free energy penalty when the phases mix 
and the cluster surfaces become fluffy. This is quite sim- 
ilar to the distinct behavior of cluster surfaces under nu- 
cleation versus spinodal decomposition. The suggestion 
from figures El and El is then that the dynamics changes 
the transition scenario to spinodal for all our q-values. 
In these cases the heights of the peaks are quite similar 
to those which we find for the equilibrium peaks of the 
q = 2 and q = 4 second order transitions. 



E 



9 

8 

7 

6 

5 

4 

3 

2 

1 


0.6 



80-1 
80-2 
80-4 
80-8 
80-16 




0.7 0.8 



0.9 

P 



1.1 



1.2 




FIG. 7: The largest cluster surface for the 10-state Potts 
model on 80 x 80 lattices for the rip values indicated by the 
extensions to the lattice size. 



FIG. 8: The largest cluster surface from np = 1 dynamical 
simulations on 100 x 100 lattices. Compare with the equilib- 
rium results of figure E] 



lattice. For increasing np we observe a slight decrease of 
the (3 — > /3 m j n peaks, while the f3 — > /3 max peaks increase. 
Although the peaks of the cooling and heating half-cycles 
approach one another in this way, each process is still far 
away from equilibrium as a comparison with the height 
of the q = 10 equilibrium peak of figure El shows. The 
long tails of the peaks of the (5 — > /3 ma x half-cycle de- 
crease rather rapidly with increasing np, so that S max (f3) 
approaches its equilibrium value for f3 > /Si^. 

For q = 2 (no figure shown) an approach of both peaks 
to the equilibrium peak of SVax is observed, whose height 
is for q = 2 only about 10% smaller than the height of the 
np = 1 dynamical peak. We take this as an indication 
that in the range of our dynamical speeds the phase con- 
version mechanism is always spinodal, independently of 
the order of the equilibrium transition. Figure El makes 
this point by contrasting the equilibrium results of fig- 
ure El with the n/3 = l dynamical results. In the next 
subsection we analyze our structure functions data with 
respect to this scenario. 

The locations of the equilibrium peaks are closer to 
the /3~ cak values of the dynamical /3 — > /3 m in heating half- 
cycles than to the (3p Cak values of the dynamical (3 — > 
Anax cooling half-cycles. This is particularly clear for 
q > 4. Our understanding of this is that the relaxation is 
faster for the heating than for the cooling half-cycle. This 
observation goes hand in hand with the interpretation of 
the higher peaks in figures El 13 an d El as being due to 
faster response times of the systems. 



Structure Functions 



The question emerges, how fast is the equilibrium sce- 
nario approached when the speed of the dynamics slows 
down? In figure [7| we plot for np — 1, 2, 4, 8, and 16 
our S max (f3) results of the 10-state model on an 80 x 80 



During our simulations we recorded the structure func- 
tion for the following momenta: 

ki = {2nL- 1 ,0) and (O^ttZT 1 ) (14) 
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k<2 



(2ttL^ 1 ,27tL- 1 ) 
(47rL -1 ,0) and (0,47rL -1 ) 
{AttL- 1 ^ttL- 1 ) and (27TL -1 , iwL' 



(15) 
(16) 
- 1 ) (17) 
(18) 



The structure functions are averaged over rotationally 
equivalent momenta. In the following we use the notation 
Ski — Ski(fi), {i = 1,...,5) for the structure function 
S(k, t), when the vector k is fcj and the time dependence 
is dictated by (3 — [3{t). Spinodal decomposition is char- 
acterized by an explosive growth in the low momentum 
modes, while the high momentum modes relax to their 
equilibrium values. 
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FIG. 9: The structure function S kl {0) for the 10-state Potts 
model and np = 1 dynamics. 
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FIG. 10: The structure function S kl (P) for tne 2-state Potts 
model and np = 1 dynamics. 



FIG. 11: The structure function S kl (/3) from np = 1 dy- 
namical simulations on 80 x 80 lattices together with their 
equilibrium values. 
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FIG. 12: The structure function S kl for the 10-state Potts 
model on 80 x 80 lattices for dynamical simulations with np 
as indicated by the extensions to the lattice size. The inlay 
of the left enlarges the equilibrium peak together with the 
heating data. For the inlay /3 is mapped on 0.5 + 2 * (f3 — 0.66) 
and S kl on 0.002 + 4 * S kl ■ 



In figures I§1 and 1101 we give the q = 10 and q = 2 re- 
sults for Skj^iP)- The hysteresis flow is indicated by the 
arrows. Figure ITTI shows the g-dependence of the n« = 1 
dynamics together with our equilibrium results. In com- 
parison with the S max (0) figures 03 H3 and |H1 several 
differences and similarities deserve to be mentioned: 

1. For our np values the Sk t structure functions have 
a pronounced maximum on the (3 — > /3 ma x half- 
cycle. This is cooling in the spin system language 
used in this paper and heating (i.e., confinement to 
deconfinement) in an analogue QCD system. 

2. Different ordinate scales are chosen in figures 
and because the magnitudes of the peaks show 
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so high that the peaks under heating become overshad- 
owed. This can be made explicit by first equilibrating 
the systems at /3 max . Then peaks of similar size as the 
equilibrium peaks appear on the /3 max — * Pmin half-cycle. 
For our slower dynamics the (/?) structure function 
peaks for the (3 — ► f3 m - ln half-cycle become visible without 
equilibrating first at /3 max . For increasing np they ap- 
proach the equilibrium peaks (see the inlays of figures IT21 
and 1 1 3|) . For both half-cycles the approach to equilib- 
rium appears to happen only for a really slow dynamics. 
Very CPU time consuming simulations of np 3> 16 values 
would be needed to follow this in detail. 



2. ki, i > 1 and quenching 



FIG. 13: The structure function Ski for the 2-state Potts 
model on 80 x 80 lattices for dynamical simulations with np 
as indicated by the extensions to the lattice size. The inlay 
of the left enlarges the equilibrium peak together with the 
heating data. For the inlay /3 is shifted by —0.13 and Sk x by 
+0.002. 



a considerable g-dependence. That is exhibited in 
figure ITT1 

3. As in figure [5] for S max the equilibrium peaks of 
Sk ± increase from q — 10 to q — 2 (see figure ITT|l. 
However, increasing np from 1 to 16, the approach 
of the Sfcj (13) function to their equilibrium is even 
for j3 > @y Ca k rather slow. This is shown for q = 10 
in figure 1121 To a large extent it holds still for 
g = 2, where for [3 > /3p eak the equilibrium appears 
to become approached for np = 16, as is shown in 
figure ED 

4. The inlays of figures IT51 and IT31 enlarge the equilib- 
rium peaks together with the heating (/? — > /3 m ; n ) 
data. For q = 10 as well as for q — 2 we find that 
the heating data develop a peak with increasing 
np, which may eventually approach the equilibrium 
peak. However, in both cases it appears that the 
heating peak wants first to merge with the cooling 
peak. For q — 10 the np — 16 cooling peak is much 
larger than the equilibrium peak and the heating 
peaks are all the time increasing. Eventually, both 
peaks should start to decrease towards the equilib- 
rium data. For q = 2 the cooling peaks decrease 
rapidly with increasing np and the np — 16 cool- 
ing peak undershoots the equilibrium peak. Cool- 
ing and heating peak move towards merging and 
should approach the equilibrium peak from below. 

The Sfej structure function peaks strongly under cool- 
ing and less under heating. This is presumably related to 
the fact that the spin variables get ordered at low tem- 
peratures (like the Polyakov loops get ordered at high 
temperatures), thus allowing for g-ality order-order do- 
mains. For our fast dynamics the Sk t values at /3 max are 




FIG. 14: Hysteresis of the Sk(t) structure functions for the 
10-state Potts model on an 80 x 80 lattice (np — 1). 
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FIG. 15: Hysteresis of the Sk(t) structure functions for the 
2-state Potts model on an 80 x 80 lattice (np = 1). 

Miller and Ogilvie 4] investigated the dynamics of 
SU(2) gauge theory after quenching from a low to a high 
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FIG. 16: Time evolution of the Sk(t) structure functions for 
the 10-state Potts model on an 80 x 80 lattice after quenching 
from m in = 0.4 to /3 max =0.8. 
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FIG. 17: Time evolution of the S k (t) structure functions for 
the 2-state Potts model on an 80 x 80 lattice after quenching 
from /3 min = 0.2 to /3 max = 0.6. 

physical temperature (corresponding to the /3 m in 
half-cycle of the spin system). They report a critical 
value k c , so that modes grow (do not grow) exponentially 
for k < k c (k > k c ). 

In figures n~T1 and IT51 we show for our rip = 1 dynamics 
all structure functions, which we have measured. It is 
notable that we observe a large gap between the peaks 
for S ki (l3), i = 1,2 and for S k] {[3), j > 3. We also per- 
formed quenching runs. In figure 1161 we show the time 
evolution on an 80 x 80 lattice after quenching the 10- 
state Potts model from (3 — 0.4 to (3 = 0.8 and in figurclTTI 
we show the time evolution after quenching the 2-state 
Potts model from (3 = 0.2 to (3 — 0.6. As in our hystere- 
sis investigations averages of 640 independent repetitions 
are taken. Again, we find a large gap between the peaks 
for S ki {P), i = 1,2 and for S kj ((3), j > 3. Further it is 
remarkable that the height of the peaks in the hysteresis 
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FIG. 18: Time evolution of the Sk(t) structure functions for 
the 2-state Potts model on an 80 x 80 lattice after quenching 
from /3 m in = 0.2 to /3 max = 0.6 at a magnetic field h = 0.01. 



cooling half-cycles and under quenching are almost iden- 
tical, while the time scale is according to equation Q 
extended to 3200 sweeps for the hysteresis curve. 

For figure ITTjl we have changed the second order phase 
transition of the 2-state Potts model to a crossover by 
adding the term hJ2f^a-(r,t),q with a small magnetic 
field h, h = 0.01, with respect to the spin qo in the Hamil- 
tonian Q . In essence the time evolution after quenching 
is similar as without the magnetic field, only that the 
magnitude of the S kl and S k2 peaks decreases, maintain- 
ing still a clearly visible gap to the peaks of the S kj , 
j > 3 structure functions. Increasing the magnetic field 
further, to h = 0.1, the large peaks disappear altogether 
by merging into the small peaks. The signal of a transi- 
tion is possibly lost for such high values of the magnetic 
field. 

The observed peaks suggest a critical fc.; value between 
i — 2 and i = 3 for the Potts models. However, the issue 
is more subtle. Cahn-Hilliard theory 26, 40] (model B for 
reviews see |2(], |2^, |2£| ) predicts an exponential growth 
of the low momentum structure function in the initial 
part of the time evolution after quenching (also applies to 
model A). Whether such an exponential growth is found 
or not was used by Miller and Ogilvie |4| to determine 
the critical k c between the low and the high momentum 
mode. However, none of our structure functions in fig- 
ures El to El shows initial exponential growth. This is 
already kind of obvious by looking at the figures, where 
the shape of the increasing parts of the curves is always 
concave and is quantitatively determined by performing 
fits. A likely explanation is that the Cahn-Hilliard the- 
ory relies on approximations which are not justified in 
our 2d models. In 3d we find exponential growth in the 
very early stage of the time evolution after quenching 
(incrementing then (3 after every update and fitting the 
time evolution within the first sweep [13|). Interestingly, 
our hysteresis curves of figures El an d 1151 which rely on 
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the smoother dynamics of temperature changes in small 
steps, show an initial exponential growth. 

We find no peaks in the structur factors when we 
quench from an ordered initial state into the disordered 
phase. The reason may be that the structure factor is 
defined with respect to the order parameter. Under a 
quench into the ordered phase g-ality order-order domain 
may emerge, whereas there is only one disordered phase. 

IV. SUMMARY AND CONCLUSIONS 

Our energy hysteresis method allows for dynamical es- 
timates of the equilibrium transition temperatures for 
first as well as for second order phase transitions. While 
the precision of these estimates is not competitive with 
those of equilibrium investigations, the hysteresis method 
provides information about dynamically rooted devia- 
tions of accompanying physical observables from their 
equilibrium values. For second order transitions we find 
that the dynamics generates a latent heat and for a weak 
first order transition we find a 'dynamical' latent heat 
much larger than its equilibrium value, whereas for a 
strong first order transition the dynamical latent heat 
agrees with the equilibrium value (the magnetization al- 
lows for a similar analysis, which is not reported here). 

In our analysis of 2d Potts models, we find spinodal 
decomposition to be the dominant feature as soon as 
we turn on the dynamics. For instance, the equilibrium 
(quasi static) phase conversion of first order phase tran- 
sitions is due to nucleation. Even our slowest dynamics 
(np = 16) changes the phase conversion of the investi- 
gated weak (q = 5) and strong (q = 10) first order tran- 
sitions from nucleation to spinodal. For the q = 2 (weak) 
second order transition the np = 16 dynamics appears to 
be already rather close to equilibrium, which is formally 
reached for np — > oo. These results are mainly based 
on analyzing the dynamical time evolution of Fortuin- 
Kasteleyn (FK) clusters and structure functions. 

• For FK clusters we find that the largest cluster sur- 
face area is quite sensitive to dynamical effects and 
yields for all considered q- values signals in favor of a 
spinodal decomposition on the cooling (f3 — > /3 max ) 
and heating (f3 — > /3 m j n ) half-cycles of our hysteresis 
loops. This may be illustrated by comparing the re- 
sults of our fast (np = 1) dynamics of figure |S] with 
the equilibrium results of figure [BJ For the first 
order transitions the dynamics enhances the peak 
values to take on similar values as one finds for the 
q = 2 and q = 4 second order configurations near 
the critical point. 

• For the structure factor our np = 1 dynamics leads 
on the /3 m i n — > P max half-cycle to amplitude max- 



ima, which are considerably larger than those from 
the second order equilibrium configurations, see fig- 
ure The dynamical peaks on the /3 max — ► Pmin 
half-cycle are comparable to those of from the equi- 
librium configurations, which has for figure ITT1 the 
consequence that they are not visible at all, because 
the systems are still out of equilibrium at /3 max and 
on their return path. That is of potential interest 
for heavy ion collision, where our /3 m i n — > /3 max half 
cycle corresponds to heating for which the dynam- 
ics of the experiment is definitely fast. We have 
no entirely satisfactory theoretical explanation for 
the observed asymmetry, but think that is is re- 
lated to the fact that g-ality order-order domain 
may emerge in the ordered (deconfined) phase. 

Using quenching, we find dynamical signals surviving 
even after the proper (second order) phase transition is 
converted into a crossover. Moving then far away from 
the transition, the dynamical signals fade away too and 
the issue of crossovers requires further investigations. 

Our computer programs allow to extend the present 
study to the 3-state Potts model in three dimensions 
with an external magnetic field representing quark ef- 
fects. In a more remote future, one could carry out simi- 
lar studies for quenched and even full QCD. But none of 
these studies could resolve the problem of a quantitative 
relationship between the Glauber time scale of our Eu- 
clidean dynamics and the time scale of the Minkowskian 
dynamics in the real world. In this context it is of in- 
terest that Pisarski and Dumitru developed recently a 
Polyakov loop model 0] which allows for simulations in 
the Minkowskian formulation 0]. It mav be possible 
to address questions similar to those raised in our paper 
within the hyperbolic dynamics of their model. 

Most likely the aim of such studies cannot be to make 
precise quantitative predictions. Instead, one may have 
to be content with illustrating the effect of different 
speeds of the phase conversion on the observable sig- 
nals qualitatively. If spinodal decomposition of Polyakov 
loops is indeed realized in heavy ion collisions, one may 
observe an enhancement in the production of low-energy 
gluons. 
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